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COMPUTATIONAL  ASPECTS  OF  CONSTRAINED  ESTIMATION 


INTRODUCTION 


The  representation  of  physical  systems  by  logical -mathematical  models  and 
the  subsequent  realization  of  these  models  as  computational  programs  dominates 
engineering  and  scientific  activities.  Generally,  a  model  characterizes  a 
cause-effect  relationship  as  a  mapping  of  a  set  C  into  a  set  E  flm:C  >  E).* 
Problems  concerned  with  such  mappings  are  variously  classified  depending  upon 
the  information  available  and  the  information  sought.  The  direct  (normal  or 
analysis)  problem  seeks  to  generate  the  set  of  effects  from  a  set  of  causes, 
while  those  problems  concerned  with  reversing  the  cause-effect  relationship 
are  termed  inverse  or  indirect  problems.  Synthesis  (identification)  problems 
require  the  determination  of  the  laws  or  mappings  that  govern  cause-effect 
relationships.  A  fundamental  concern  is  the  awareness  of  conditions  for 
which  solutions  to  such  problems  are  revealed. 

A  problem  is  said  to  be  "well-posed"  in  the  sense  of  Hadamard*  if  a  unique 
solution  exists  and  arbitrary  small  changes  in  model  parameters  lead  to 
correspondingly  small  variations  in  the  solution. This  latter  condition 
of  stability  is  significant  considering  that  it  is  quite  standard  to  assure 
the  existence  and  uniqueness  of  solutions  and  yet  accommodate  the  possibility 
of  results  having  little  physical  relevance.  The  viewpoint  is  maintained  that 
physical  systems  that  lead  to  problem  statements  having  potentially  unstable 
solutions  are  not  unusual  and  that  methods  for  the  appraisal  and  resolution  of 
the  problem  exist  and  add  little  to  the  computational  burden  of  the  problem. 


*A  further  characterization  of  Hadamard's  conditions  is  given  by 
27 

Rutman  where  for  A,  a  linear  map  in  Banach  space,  b  €  ImA  (existence)  and 
kerA  »  0  (uniqueness)  imply  an  algebraic  well-posedness,  and  ImA  *  ImA 
(continuous  dependence  of  the  solution)  implies  a  topological  well-posedness. 


1 


Inverse  problems  comprise  the  majority  of  those  problems  that  may  be 

considered  ill-posed,  viz.,  unique  solutions  that  are  unstable.  Usually,  the 

analysis  of  direct  problems  of  quadrature  and  problem?  of  synthesis  are  easily 

reformulated  as  inverse  problems.  Examination  of  the  causes  of  i 1 1-posedness 

suggests  the  following  classifications:  problems  that  are  shown  to  be 

3 

ill-posed  as  stated  (e.g.,  by  arguments  of  the  Riemann-Lebesque  lemma  ); 

5 

problems  having  parameters  contaminated  by  noise;  and  problems  for  which 
computational  solutions  are  sought  and  thus  vitiated  by  errors  of  truncation. 
Most  often,  il 1-posedness  results  from  combinations  of  these  factors,  which  in 
some  cases  provide  illusions  of  reliable  analysis,  but  in  all  cases  compound 
initial  difficulties.6 

Instances  of  ill-posed  problems  are  frequently  encountered  in  applied 
situations,  numerous  examples  of  which  can  be  found.1 ,3’4’7-10  Often, 
these  applications  require  solving  systems  of  simultaneous  linear  and 
non-linear  equations,  ordinary  and  partial  differential  equations,  and 
integral  equations.  Of  contemporary  interest  are  industrial  applications 

that  include  inverse  scattering  problems,11-14  ocean  acoustic 

IB  16  17  23 

tomography,  seismology,  inverse  wave  propagation,  image 

restoration  (typically  in  problems  requiring  restoration  beyond  the 

24  25  26 

diffraction  limit),  mathematical  optimization,  data  query,  and 

25  27—30 

optimal  control.  *  A  common  subproblem  of  continuing  mathematical 

interest  is  the  solution  of  the  general  linear  problem,  which  is  seen  to  be  a 

finite-dimensional  representation  of  the  Fredholm  integral  equation  of  the 

31 

first  kind.  In  the  case  of  non-linear  problems,  linearization  techniques 

32 

or  other  innovations  are  often  successful  in  reducing  the  problem  to  a 
linear  equivalent.  A  main  consideration  is  the  investigation  of  computational 
methods  appropriate  for  the  solution  of  linear  ill-posed  problems  and  the 
application  of  methods  of  constrained  estimation. 


2 


GENERAL  LINEAR  PROBLEM 


Problems  for  the  solution  of  Fredholm  integral  equations  of  the  first  kind, 
that  is. 


n2 

f  K(t,$)f(s)ds  =  g(t)(c1  <  t  <  c2)  (1) 

occur  frequently  whenever  input  data  f(s)  are  to  be  determined  from 

measurements  g(t)  obtained  from  the  output  of  a  device  or  instrument;  the 

kernel  function  K(t,s)  categorizes  the  measurement  process.  As  a  practical 

example,  consider  the  problem  of  estimating  the  acoustic  field  N(e),  from  beam 

33 

measurements  M(y)  of  a  line  array  consisting  of  k  equally  spaced  elements. 

This  problem  results  in  the  expression 

sin  jkirc /t [cos (©  i>)  -  cos(Y)]}  N(e)d6! 

ksin2  {itC/t[cos(©  -<!>)-  cos(y  m 

0 

where  ©  is  the  azimuth  measurement,  y  corresponds  to  the  beam  steering  direc¬ 
tion,  i k  is  the  array  heading,  t  is  the  wavelength  of  the  input,  and  c  is  the 
array  element  spacing.  For  continuous  kernels,  this  problem  is  ill-posed.*® 

Algebraization  of  (1)  proceeds  by  application  of  an  appropriate  rule  of 
quadrature  w. ,  which  results  in  the  approximation 

J 

m 

]T  =  wjK(ti,s.)f(sj)  =  g(t1)  (i  -  1,2, — ,n)  . 
j  *  0 

A  typical  selection  of  abscissae  Sj  is  given  by  Sj  =  a  +  h(j  -  1),  where 
the  displacement  h  =  (n^  -  n^)/(m  -  1).  For  a  Simpson's  rule  of 
quadrature, 


3 


h/3  (for  j  »  1  or  m) 

[3  +  (-l)^]h/3  (otherwise) 

and  m  is  odd.  There  then  remains  the  selection  of  the  quantization 

parameters  n  and  m,  which  are  chosen  to  be  small  enough  to  reduce  the  degree 

of  numerical  truncation,  but  large  enough  to  assure  adequate  resolution  and 

34 

representation  of  the  physical  process.  In  Hunt  this  dilemma,  encountered 
in  the  naive  solution  of  equation  (1),  is  demonstrated  and  further 
characterizes  the  il  1-posedness-  of  such  problems. 

A  compact  representation  of  a  finite-dimensional  linear  system  is  provided 
by  a  matrix  expression,  which  for  (1)  results  in 

Ax  =  b  .  (2) 

Conditions  for  the  existence  of  solutions  to  (1)  follow  by  Picard's 
criteria.  Alternatively,  equation  (2)  admits  a  solution  if  and  only  if 
the  coefficient  matrix  [A]  and  the  augmented  matrix  [A : b]  are  of  equal  rank; 
i.e.,  the  vector  b  is  a  linear  combination  of  the  columns  of  A,  in  which  case 
the  linear  system  is  consistent.  If,  in  addition,  the  matrix  A  is  square  and 
of  full  rank,  then  the  solution  to  (2)  is  unique  and  can  be  generated  by 
application  of  Cramer's  rule  (or  the  inversion  of  A).  For  an  over determined 
linear  system,  unique  "solutions"  are  provided  by  transforming  equation  (2) 
into  a  system  of  normal  equations 

A' Ax  =  A'b  ,  (3) 

which  is  seen  to  be  mathematically  equivalent  to  the  linear  least-squares 
problem. 

But  Cramer's  rule  is  applicable  only  if  the  determinant  det(A)  of  the 

square  coefficient  matrix  A,  in  (2),  is  non-zero.  In  situations  where  the 

determinant  is  zero,  the  matrix  is  said  to  be  non-regular  or  singular.  In 

linear  systems  having  small  determinants,  A-*  contains  elements  of  large 

36 

amplitude  that  amplify  small  variations  in  the  measurement  vector  b. 
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This  deterioration  in  the  conditioning  of  the  linear  system  (ill-conditioning) 

is  seen  as  the  numerical  manifestation  of  i 11-posedness  and  is  manifested  by 

increasing  interdependence  among  the  rows  of  the  coefficient  matrix  or 

smoothness  of  the  kernel  in  equation  (1).  Symmetrization  only  exacerbates 

2  6 

this  condition  since,  for  A  (a  square  matrix),  det(A'A)  =  det(A)  . 

Adequate  appraisal  of  conditioning  is  then  seen  as  prerequisite  to  the 
generation  of  solutions  to  linear  systems. 

A  useful  method  of  gaging  the  sensitivity  of  a  linear  system  is  provided 
by  the  relative  condition  number  K  where 

K  =  II  A||  ||  A_1||  .  (4) 

Typically,  the  norm  of  the  real  square  matrix  A  corresponds  to  the  maximum  of 

the  sums  of  the  moduli  in  the  rows  of  A.  The  matrix  A  represents  a  linear 

mapping,  or  transformation,  of  an  arbitrary  vector  x  into  the  vector  Ax. 

The  norm  of  A  is  then  a  measure  of  the  distortion  under  the  linear 
37 

transformation.  By  a  simple  perturbation  analysis  of  equation  (2), 
ll<sx(|  <  K  ||6b[| 

lixir  “  uw 

by  which  for  large  K  it  is  seen  that  small  relative  variations  in  A  and  b  are 
magnified  in  x,  a  restatement  of  the  condition  of  i 11-posedness.  However, 
determining  equation  (4)  is  complicated  by  the  uncertainty  in  the  estimate  of 
||A-*II  when  K  is  large. 

Alternatively,  the  condition  number  can  be  determined  by  the  square  root 

op 

of  the  ratio  of  the  largest  eigenvalue  of  AA'  to  its  smallest  eigenvalue. 

It  is  then  apparent  that  ill-conditioning  is  associated  with. eigenvalues  close 
to  zero.  A  dominant  eigenvalue  x  of  D  «  AA',  where  D  is  similar  to  a 
diagonal  matrix,  can  be  determined  from  the  sequence 

w  =  Dw_  ,/xm  , 
m  m-1  m  * 

where  w„  is  an  arbitrary  vector  and  x„  is  the  maximum  element  of  Dw 

o  m  m-l 


Similarly,  the  sequence 


U=  (0  -  XmI)un  i / X 
n  m  '  n-l  n 

provides  an  estimate  of  the  minimal  eigenvalue  xn*  For  non-diagonal  type 

matrices,  deflation  procedures  may  be  appropriate  for  determining  the  minimal 
39 

eigenv  alue. 
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GENERALIZED  INVERSE  SOLUTION 


Often,  the  occurrence  of  ill -posed  problems  in  practice  motivates  the 

application  of  multiple  precision  computations  generating  costly  failures 

and  a  skepticism  in  the  promise  of  computing  machinery.  Alternative  remedies 

(e.g.,  preconditioning  via  scaling^),  the  use  of  orthogonal  (or  nearly 

orthoqonal )  base  functions, ^  large  number  arithmetic,  interval 
42-45  46  47 

arithmetic,  direct  search  methods,  ’  and  forward/ backward  error 
48-50 

analysis  may  provide  some  insight  into  the  characteristics  of  the 
solution,  but  usually  at  excessive  cost.  Even  problems  of  moderate 
dimensionality  can  rapidly  exhaust  available  memory  capacities  as  machine 
truncation  errors  are  compounded.  Invariably,  information  is  lost  as  real 
numbers  (viz.,  irrational  numbers)  are  mapped  onto  the  sieve-like  range  space 
that  persists  in  the  computational  environment.  Similar  effects  result  as 
measurements  become  contaminated  by  noise.  Attempts  to  ameliorate  the  effects 
of  lost  information  are  frustrated  unless  additional  information  is  furnished. 

Although  existence  of  a  solution  to  equation  (3)  is  assured  for  an 
arbitrary  coefficient  matrix,  uniqueness  is  provided  if  and  only  if  the 
coefficient  matrix  is  of  full  rank.  Since  the  set  of  all  least-squares 

3 

solutions  forms  a  closed  convex  set,  a  unique  element  can  be  selected  that 

solves  equation  (3).  An  appropriate  selection  is  given  by  the  least-squares 

solution  of  equation  (2),  which  is  of  minimum  norm  variously  denoted  by  the 

51 

generalized  inverse  or  pseudo-inverse  solution.' 

For  an  arbitrary  m  x  n  matrix  A,  the  generalized  inverse  is  defined  by  the 
unique  n  x  m  matrix  A  ,  which  satisfies  Penrose's  lemmas: 

A+AA+  =>  A+ 

AA+A  =  A 
(AA+)'  -  AA+ 

(A+A)'  -  A+A  . 
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Computationally  advantageous  is  the  representation  of  A+  in  terms  of  the 
singular-value  decomposition  of  A,  that  is, 

A  =  UL1/2V'  .  (5) 

The  matrices  U  and  V*  result  from  columns  formed  by  the  eigenvectors  of  AA' 

and  A'A  and  L  is  a  m  x  n  matrix  composed  of  a  k  x  k  (k  =  rank  (A))  diagonal 

12 

matrix  of  the  corresponding  eigenvalues  with  the  remainder  zero-filled. 

The  generalized  inverse  is  then  given  by 

A+  =  V 1 L-1/2U '  .  (6) 

By  Lagrange  minimization,  the  least- squares  solution  of  (2)  having  minimum 
norm  results  in  the  expression 

A+u  =  (A'A  +  U2I )-1A'  .  (7) 

Substituting  (5)  into  (7),  and  noting  that  U  and  V  are  unitary,  results  in 
A+p  =  V'L1/2[L  +  y2I ]-1U'  , 

which  is  equivalent  to  equation  (6)  as  y  approaches  zero. 

Equation  (7)  corresponds  to  the  ridge  inverse  (compare  the  Levenberg- 

CO  CO 

Marquardt  procedure  and  damped  least-squares  )  as  applied  in  ridge 
54 

regression,  or  the  constrained  estimate  of  the  solution  of  (2)  and  is  seen 

55 

to  be  equivalent  to  an  approximate  generalized  inverse.  The  ridge 

estimate  may  be  seen  to  be  a  type  of  weighted  average  between  the  input  data 

56 

and  supplemental  information.  This  idea  can  be  further  extended  to 
include  varied  qualities  of  information  that  conform  to  the  anticipated  nature 
of  the  solution. 
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GENERALIZED  CONSTRAINED  ESTIMATION 


The  application  of  supplemental  information  to  the  solution  of  ill-posed 
problems  is  quite  common.  For  example,  Wiener  theory  requires  statistics  con¬ 
cerning  the  signal  and  noise  processes. ^  Such  quantitative  inputs 
(a  priori  information),  when  available,  may  not  be  sufficient  for  the  solution 
of  ill-posed  problems;  additional  information  will  often  be  required. 

Generally,  all  relevant  factors  characterizing  the  nature  of  the  solution 

27 

may  be  described  as  exogenous  information.  Because  of  the  variety  of 
information  that  may  influence  the  problem,  generalized  solution  methods  are 

CO 

rare.  An  early  approach  by  Kreisel  seeks  smooth  solutions  to  equation  (1) 
utilizing  exponential  weighting  functions  that  effectively  neutralize  the 
severe  oscillations  prevalent  in  ill-posed  solutions.  Moreover,  smooth 
solutions  are  readily  justified  as  representative  of  most  physical  processes. 

cq  cn 

Similar  constraints  motivated  the  development  by  Philips  and  Twomey 
for  equation  (1)  in  the  form 

min  { || Ax  -  bll  +  a<t>2(x) }  ,  (8) 

where  <t>  (x)  corresponds  to  (side)  conditions  imposed  on  the  least-squares 
solution,  and  a  >  0  is  the  degree  to  which  such  conditions  should  influence 
the  solution.  For  smoothness  constraints,  <j>  (x)  *  II  »  where  4>  corre¬ 
sponds  to  a  second-difference  approximation  to  the  solution  and  is  of  the  form 


9 


A  general  solution  to  (8)  is  given  by 


X  =  ( A'A  +  ail)  1  ( A ' b  +  aC)  , 

where  c  corresponds  to  a  known  (absolute)  bias  and  a  =  4>'$.  The  potential  of 
using  alternative  constraints  was  also  suggested  by  Twomey  (e.g.,  (x.ftx)  is 
equivalent  to  the  variance  of  the  solution  where 


(n-l)/n. 

i  =  j 

-1/n, 

i  *  J 

(i,j  -  1,2, ...,n) 

and  where  n  is  the  number  of  measurements.  Constraints  in  terms  of  entropy 
functions  can  be  found  in  Smith. 


REGULARIZATION 

The  approach  indicated  in  (8)  was  also  independently  investigated 
63 

by  Tikhonov  and  resulted  in  his  formalizing  the  method  of 
regularization.*  By  regularization  is  meant  an  adjustment  of  the  initial 
problem  that  admits  proper  solutions.  Justification  of  the  method  is  given 
in  various  sources^’^’^-^  where  a  general  form  of  the  smoothing 


functional  in  (8)  is  given  by 

1  P 

'  dJ‘x(e)  " 

E  ¥■> 

0  j  =  0 

.  de^  . 

Thus,  for  aQ  =0,  a^  =  0,  a?  =  1,  and  p  =  2,  there  results  the  smoothing 
functional  corresponding  to  (9)  and  in  a  similar  manner  the  estimation  of 
equation  (7)  is  defined  for  aQ  =  pQ  =  0. 


*Also  called  the  Tikhonov-Mil ler  method. 


10 


The  method  of  regularization  is  given  a  number  of  statistical 
interpretations^’67-77  in  which  the  Tikhonov  method  is  supplemented  by  a 
Wiener  technique  or  a  Bayesian  strategy.  An  estimate  for  the  problem  of 
equation  (9)  can  be  given  by 

X  -  (A'WA  +  aH)-1A'Wb  ,  (10) 

Qt 

where  W  is  a  weighting  coefficient  matrix  generated  from  known  statistics 

23  78 

concerning  the  measurements  b.  In  Edenhofer  and  Varah  H  corresponds 
to  the  assumed  known  covariance  matrix  of  the  data  x,  and  W  corresponds  to  the 
covariance  matrix  of  the  measurements  (a  particular  form  of  the  weighting 
matrix) . 

A  problem  that  remains  is  the  selection  of  the  regularization  parameter  in 
(8).  This  formulation  is  seen  to  correspond  to  the  Lagrange  minimization  for 
the  problem: 

m1n{*(xa):  ||Axa  -  b  II  2  =  ||  e  ||  2  j  ,  (11) 

*a 

where  the  vector  e  represents  the  measurement  error.  A  suitable  a  is  then 

chosen  so  that  the  equality  in  (11)  is  satisfied.  Iterative  procedures  for 

57  74  75 

determining  a  then  require  an  estimate  of  e.  *  ’  To  facilitate  such 

iterations,  the  estimates  of  equation  (10)  may  be  obtained  by 

xfl  -  [I  +  aH(A’WA)-1]-1x0  ,  (12) 

which  may  be  of  computational  advantage. 

78 

Varah  suggests  an  interactive  graphics  approach  for  the  solution  of 
the  basic  problem  of  equation  (1).  Graphical  display  of  a  problem  may  be 
successful  because  it  motivates  the  analyst  to  employ  perceptual  abilities 
that  the  alphanumeric  display  precludes.  This  approach  can  be  extended  to 
provide  a  subjective  method  for  approximating  a,  as  well  as  applying 
additional  regularization  criteria. 
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Oirect  methods  for  estimating  o  are  often  desirable;  but  without 

12 

adequate  noise  estimates  they  are  less  accurate.  Labianca  provides  an 
empirical  approach  where  a  is  selected  as  a  fraction  of  the  maximum 
singular-value  of  A'A,  the  fraction  0.0136  being  a  suitable  choice. 

QUADRATIC  OPTIMIZATION  METHODS 

For  many  physical  problems  a  non-negative  solution  (e.g.,  probability 
distributions,  signal  spectra)  is  required.  An  n-step  method  for  selecting  a 
regularized  solution  given  non-negativity  constraints  is  offered  by 
Turchin.7^  This  basic  approach  is  further  exploited  by  Rutman®^  to 
include  additional  qualities  of  constraints  in  a  selective  manner  if  required. 

Solution  of  the  general  linear  problem  (2)  with  non- negativity  constraints 
is  equivalent  to  the  problem 

minjz(x)  =  l/2(x,Dx)  -  (q,x):x  >  oj*  ,  (13) 

where  D  =  A'A  and  q  *  A'b.  A  regularized  solution  requires  that  D  *  A'WA  +  oH 
and  q  =  A'Wb.  For  the  situation  where  the  solution  x°  to  the  unconstrained 
problem  is  negative,  a  proximal  non-negative  solution  xm  to  x°  exists** 
and  is  determined  by  the  following  algorithm: 


♦Equivalently,  Z(x)  =.  l/2(Dx  -  2q,x). 

**A  solution  at  the  origin  is  always  possible. 


Algorithm  I 


Step  0:  Given  D  and  x°,  then  q  *  Dx°,  S  =  |i:x°  >  0}  and 
xP  =  max[0,x°],  Vi . 


r  T-1 


Step  1:  Set  x  =  T”  z,  where  T  *  [d^]  and  2  =  D *  i,j  € 

Step  2:  If  x':  >  0,  Vi,  then  go  to  Step  4. 

Step  3:  Set  S  =  S\{j r ^  =  min[x^/(xP  -  xT )3 ,  x?  =  x1?  +  ( 

and  go  to  Step  1. 

*  r 

Step  4:  If  e.j (Tx  -  z)  >  0,  Vi,  then  go  to  Step  6. 

Step  5:  Set  S  *  SU{jV:r  .min[e.(Txr  -  z)],  and  go  to  Step  1. 

1  1  J  •  1 


SteD  6:  Set  xm  *  xr  . 


The  vector  e^  with  ith  element  unity  and  all  others  zero  is  called  the  unit 
column  matrix. 

Additional  constraints  can  be  accommodated  for  equation  (13)  by  intro¬ 
ducing  the  transformation  Rx  =  x,  where  R  is  non-singular.  Equation  (13)  then 
becomes 


min{Z{x)  =  l/2(x,Dx)  -  (q,x):x  >  Of  , 
x 


which  is  essentially  reduced  to  the  constraints  of  the  type  of  equation  (13). 
The  utility  of  the  transformation  is  demonstrated  by  examining  several  trans¬ 
formation  constraint  matrix  pairs  R  and  R~*. 


3 


Monotonicity  constraints  often  occur  in  physical  problems  (e.g., 
polynomials  with  positive  coefficients  and  cumulative  probability 
distributions).  A  monoton ical 1 y  increasing  function  x(e)  has  the  property 
that  dx(e)/de  ^  0*  Choosing 

-1  1 
-1  1 

.  -1  l 

1  _ 


corresponds  to  imposing  the  constraint  that  the  first  differences  of  the 
restored  function  x  be  non-negative.  The  matrix 


performs  the  incremental  summations  of  the  restored  differences. 


A  similar  development  leads  to  unimodality  constraints  given  by  the 
matrices 


and 
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-1  -1 
-1 


1 

1 


1  .  .  -1 

1  .  .  -1  -1 


where  n  is  odd  and  the  row  consisting  of  the  unitary  element  corresponds  to 
the  mode  position. 


Convexity  constraints  are  imposed  by  requiring  the  second  derivative  of 
the  restored  function  to  be  non-negative.  This  leads  to  the  matrices 


R 


-1 


-2 

1 


1 

-2 


1 

•  • 
1  -2 


and 


R  =  O^]  =  [r ji ]  *  C-i (n  -  j  +  1) /( n  +  1)],  i  >  j  . 

Concavity  constraints  are  developed  in  a  similar  manner. 

Non-negativity/non-positivity  constraints  are  simply  developed  by  matrices 
of  the  type  R  »  diag|...l,...-l,...j.  Further  selectivity  is  imposed  by 
limiting  the  set  of  active  constraints  to  the  intersection  of  the  set  of 
selected  coordinates  and  the  set  of  positive  constraints  S  in  Algorithm  I. 

After  determining  xm,  a  new  value  of  the  regularization  parameter  a*  is 
recomputed  to  account  for  the  adjustment  in  the  estimate  x°  due  to  non¬ 
negativity  and  additional  constraints.  In  Turchin6*  the  method  for 
recomputing  the  regularization  parameter  is  given  by 
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V 


JB*~*y*L>*J~*>  «>■ '  -■ 


a*  ■  a(n*/n)^  , 

where  n*  is  the  number  of  elements  in  the  estimate  xm  that  differ  from  zero. 
This  new  value  of  the  regularization  parameter  is  then  used  to  compute  an 
adjusted  estimate  xm  using,  for  example,  equation  (12). 

It  is  sometimes  observed  that,  although  a  finite  algorithm  is  designed  to 
terminate  with  a  solution,  in  practice  convergence  to  a  solution  does  not 
occur.  Such  is  the  case  with  Algorithm  I.  The  property  is  exhibited  that  a 
problem  may  be  mathematically  well-posed  and  remain  numerically  ill- 
conditioned.  Algorithm  I  will  require  inversion  of  matrices  of  the  order  up 
to  n-1  for  which  truncation  errors  are  prevalent.  Because  this  pathology  is 
unavoidable,  it  is  always  desirable  to  minimize  its  effects. 

The  following  iterative  algorithm  is  intended  to  solve  the  problem  of 
equation  (13)  while  avoiding  matrix  inversion: 

Algorithm  II 


Step  0:  Given  D  and  x°,  then  q  *  Dx°.  Set  i  -  0  and 
x?  =  0,  Vi  *  1,2,.. .,n  . 

Step  1:  Set  k  =  |  j:  q^/d^  -  max[q./d..],  j  *  t  J  . 

Step  2:  If  qk/dfek  <_  0,  then  go  to  Step  4. 

Step  3:  Set  i  *  k;  set  x?  =  x°  +  qk/dkk  ; 

set  q.  «  q.  -  di|(qk/dkk  ,  Vi  ;  and  go  to  Step  1. 

Step  4:  Set  xm  *  x°  . 
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Since  the  Iterative  solution  begins  at  the  origin,  a  gradient  search  (Step  1) 
selects  a  coordinate  providing  a  feasible  solution.  A  corresponding 
translation  of  the  coordinate  system  (Step  3)  adjusts  the  feasible  solution  to 
coincide  with  the  origin.  All  such  displacements  that  do  not  result  in  a 
non-feasible  solution  are  then  accumulated  to  form  the  optimum  solution. 


RECURSIVE  ESTIMATES 


In  many  applications,  computational  considerations  of  memory  size  or 
processing  time  require  alternative  formalizations  of  the  processing 
algorithms  in  order  to  realize  practical  solutions  to  problems.  Recursive 
estimates  often  provide  such  efficiencies  with  little  increase  in  algorithmic 
complexity.  For  classical  least-squares  and  the  generalized  inverse,  the 
recursive  procedure  will,  with  minor  modification,  provide  for  the  deletion  of 
particular  data  points  without  recomputation  of  all  data.  In  control 
settings,  the  Kalman  filter— a  recursive  counterpart  to  the  Wiener  filter— is 
often  used.  Recursive  interpretations  for  the  methods  of  regularization  are 
similarly  desirable. 

On-line  techniques  of  regularization  appropriate  for  identification  and 

27  67 

input  signal  recovery  have  been  previously  observed  ’  and  found  to 
require  a  redefinition  of  the  method  of  regularization.  A  simple  approach 
corresponding  to  the  recursive  approximate  generalized  inverse  is  developed 
and  is  applicable  to  problems  of  the  type  found  in  Radhakr ishna.51 

An  alternative  representation  for  the  solution  to  the  problem  of  (8)  is 
given  by  the  least-squares  solution  of 

[A*  { v6«> '] ' [ x {0] '  «  [b|0] 1 , 
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assuming  no  bias.  An  additional  measurement  b^+^  requires  a  corresponding 

addition  A^n+^  *  [a^  i",an+l  n-^  t0  the  coe^‘icient  matrix.  The 
recursive  estimate  x^n+^  is  then  given  by 

X(n+ 1)  "  x(n)  +  P(n+l)A'(n+l)(b(n+l)"A(n+l)x(n)) 

P(n+1 )  3  P(n)  "  P(n)A'(n+l) 

P(n)-(0(n)  O(n))'1  • 
where  Q|n^  =  [A 1  j  ,/S $ 1  ]  . 
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SUMMARY 


Methods  of  regularization  are  shown  to  be  extensions  of  deterministic 
least-squares  estimation.  Conventional  approaches  are  limited  by  the 
availability  of  additional  information  that  must  be  introduced  into  the 
problem  to  produce  a  solution.  Similar  limitations  hold  for  methods  of 
regularization;  however,  these  methods  have  the  advantage  of  utilizing 
qualitative  inputs  and  relaxing  quantitative  information  requirements. 

More  notably,  these  methods  provide  a  resource  for  resolving  problems  that 
are  mathematically  ill-posed  and  consequently  inappropriate  for  solution  by 
conventional  methods. 

Other  methods  comparable  with  methods  of  regularization  (e.g.,  the 

augmented  Galerkin  method  and  the  singular  value  decomposition  with  truncation 
81  82 

or  damping  ’  )  are  either  less  effective  or  prove  to  be  computationally 

80 

burdensome.  Algorithm  II  and  the  transformation  methods  in  Rutman  provide 
an  indication  of  the  utility  of  the  Philips-Tikhonov-Turchin  approach  avail¬ 
able  at  minimal  expense.  Although  industrial  application  of  regularization 
has  been  generally  limited  to  image  restoration,  additional  developments  will 
broaden  the  application  to  include  control  settings  and  the  solution  of  non¬ 
linear  problems. 
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